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The sessile microbial communities known as biofilms exhibit varying architectures as environ- 
mental factors are varied, which for immersed biofilms includes the shear rate of the surrounding 
flow. Here we modify an established agent-based biofilm model to include affine flow, and employ 
it to analyse the growth of surface roughness of single-species, three-dimensional biofilms. We find 
linear growth laws for surface geometry in both horizontal and vertical directions, and measure the 
thickness of the active surface layer, which is shown to anti-correlate with roughness. Flow is shown 
to monotonically reduce surface roughness without affecting the thickness of the active layer. We 
argue that the rapid roughening is due to non-local surface interactions mediated by the nutrient 
field, which are curtailed when advection competes with diffusion. We further argue the need for 
simplified models to elucidate the underlying mechanisms coupling flow to growth. 



PACS numbers: 87.18.Fx, 87.17.Aa, 61.43.Hv 

I. INTRODUCTION 

Biofilms are surface-associated sessile microbial com- 
munities encased in a protective polymeric matrix at least 
partly of their own production [T] [5] . Part of the healthy 
human microbiome [3l H], they can also be deleterious 
when harbouring pathogenic species and protecting them 
from biocidal treatment, such as in water distribution 
systems or medical implants [51 ISj • Biofilm architectures 
take a variety of forms, including flat, rough, rippled and 
columnar, depending on both environmental {e.g. shear 
flow, nutrient supply) and intrinsic {e.g. cell motility, in- 
tracellular communication) factors [7]-[3]. Structure can 
affect function, such as the frequently-observed channels 
that are thought to permit nutrient penetration deep into 
the fllm [lOj . A deep, quantitative understanding into 
the relationship between biofilm structure and fiow would 
therefore suggest strategies for eradicating or otherwise 
modulating biofilm formation, but no theory with pre- 
dictive capability currently exists. 

The quantitative description of the growth of rough 
surfaces, both biotic and abiotic, is an established field in 
statistical physics, in particular when the surface geom- 
etry is scale-invariant or fractal [11] . Analytical and nu- 
merical treatments of model systems have demonstrated 
that their large length-scale behaviour can typically be 
grouped into a small number of so-called universality 
classes. Which class a specific system falls into depends 
on invariant intrinsic properties, such as dimension, sym- 
metries and conserved quantities, and also the nature 
of the interactions between separated surface points, i.e. 
whether such interactions are local (strictly short-ranged) 
or non-local. In the latter case, growth at one surface 
point depends (in principle) on the current geometry of 
the entire surface. Such non-locality is known to drasti- 
cally alter the fractal surface growth picture [TTJ [T2j . 

Bacterial [13l [14] and fungal [15l [16] colonies have 
been investigated within the framework of fractal surface 
growth [T7]. However, the relevance of these findings to 



biofilms, and to which (if any) universality class biofilms 
belong, remains unclear. A recent two-dimensional study 
employing a somewhat realistic model for biofilm growth 
found complex behaviour that could not be facilely in- 
terpreted using established paradigms [181. In addition, 
none of the aforementioned models incorporate flow, de- 
spite the significant effect of biofilm architecture this is 
known to have [7]. 




FIG. 1: Snapshot of system state for 7 = 0. Particle bright- 
ness is proportional to their metabolic reaction rates Vi. The 
system size is Lx ~ Ly — 40 d™^'' and the bulk nutrient con- 
centration is Co — lOA'i/2- 

In this article, we introduce an agent-based biofilm 
model in which both the nutrient field and the biofilm 
itself is coupled to the fiow, and analyse it within the 
framework of fractal surface growth. Our hybrid scheme 
extends the Individual-based Model (IbM) P^H^ by in- 
corporating adhesive links between nearby particles, re- 
placing the ad hoc 'pushing' rules. This small but far- 
reaching extension generates a mechanically consistent 
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biofilm that can react to shear stresses apphed by the 
flow. A snapshot of our model, which we dub the mechan- 
ical IbM model or m-IbM, is shown in Fig. [T] Analysis 
reveals a rapid growth of surface roughness, both paral- 
lel and normal to the direction of mean surface growth, 
that is consistent with a linear law. This is far more 
rapid than the sub-linear scaling obeyed by canonical 
models [TT], which we attribute to a non-local surface 
interaction deriving from the long-range effects of nutri- 
ent depletion. Switching on flow, we observe a similar 
growth law but with a smaller coefficient, corresponding 
to smoother biofilms. We argue this is due to the compe- 
tition between nutrient diffusion and advection, and that 
high advection modulates the non-local surface interac- 
tions resulting in a less rough film. 

This manuscript is arranged as follows. In Sec. [ll] 
we detail the modules in our model, how they are cou- 
pled, and the algorithms employed to iterate them dur- 
ing growth. In Sec. |III| we describe analysis of growing 
films in the presence of nutrient fields of varying con- 
centrations, taking care to control the finite size effects 
that are ever-present in scale- invariant systems. Starting 
without flow, we quantify the growth of surface rough- 
ness parallel and perpendicular to the mean direction of 
growth using standard metrics, in both cases finding the 
aforementioned linear growth laws. We also relate the 
depth of actively-growing particles near the surface to the 
roughness, confirming previous findings |23j . Repeating 
the analysis with flow reveals a systematic reduction in 
roughness as the flow rate increases, while not affecting 



the thickness of this active layer. In Sec. IV we attempt to 
place our findings into the broader context of fractal sur- 
face growth. Two appendices are reserved for technical 
details. In Appendix]^ we derive analytical expressions 
for the growth of a flat film, which is used to compare 
to the numerical results. Finally, in Appendix [B] we ex- 
plain how the surface heights were extracted from our 
off-lattice simulations. 
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FIG. 2: Schematic of the model. The simulation domain con- 
sists of the biofilm and the boundary layer, which lie between 
the solid base at z = and the bulk fluid a.t z = L^. Periodic 
boundary conditions are assumed in the x and y directions. 
When flow is present, it takes the form of an afRne shear 
parallel to the a;-axis with fixed rate 7. 



II. MODEL DEFINITION 

The simulation model employed here is based on the 
Individual-based Model or IbM, which is an established 
agent-based method for the mathematical modelling of 
biofilms [T9H22] . This hybrid scheme couples discrete en- 
tities representing cells or cell aggregates to one or more 
continuous scalar fields, representing soluble factors such 
as nutrients or metabolites. In our scheme, we introduce 
a single vector fleld corresponding to the fluid velocity 
that couples to both the scalar flelds and the biofilm, 
the latter through the requirement of mechanical stabil- 
ity as explained below. We also associate a mass of EPS 
(the Extracellular Polymeric Substances that make up the 
biofilm matrix [53]) with each particle, and this is used 
to determine the elastic interactions between particles. 
We first present an overview of the central variables in 
each component of the model, before describing the time 
evolution of each in detail. A summary of the physical 
parameters and variables for each module is given in Ta- 
ble H 



A. Variables and parameters 

Our model contains three components or modules, re- 
ferred to as biomass, scalars and fluid, sharing the same 
spatial domain of a rectangular box with dimensions 
{Lx,Ly,Lz)- See Fig. [2] for a schematic diagram of the 
system geometry. The solid surface to which the biofilm 
is attached corresponds to the z = plane, and the bulk 
fluid corresponds to the upper plane z = L^- Fluid flow 
(if present) is parallel to the x-axis. Periodic boundary 
conditions are assumed in the x and y directions to avoid 
introducing wall or edge effects. 

The biomass module consists of N{t) biomass parti- 
cles i — 1 . . . N(t) at time t, each with a cellular mass 
ml and an EPS mass mf (see Fig. |3][a)). The centres of 
the particles are denoted x^. Each particle is regarded 
as spherical, with a cell diameter d1 that can be related 
to the common cell density p'^ by d1 = ^Qml/irp'^. The 
EPS associated with particle i is regarded as forming a 
spherical shell of density p° extending from the cell sur- 
face. The outer diameter of this shell is denoted and is 
related to the EPS density by d^ = ^Qml/Trp" + {d'^Y. 

For this application, the scalar module consists of a 
single scalar fleld c(x) corresponding to the concentration 
of the soluble nutrient. This enters the system from the 
bulk as per the boundary condition c(z = L^) = cq (for 
simplicity, depletion and replenishment of cq with time 
is not considered). Cells reduce the local nutrient to fuel 
their increase in mass. This reaction is regarded as lo- 
calised at the centre of each particle i, with a reaction 
rate given by the commonly-employed Michaelis-Menten 
form m, 



mi[l + Ki/2/c{xi)] 



(1) 
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Metabolic activity is converted into an increase in both 
cellular and EPS masses as 9tm^ = y^|ri| and 9tm° = 

Finally, the fluid module describes the fluid velocity 
field v(x). Here only a simple afhne shear flow is consid- 
ered, i.e. v(a;, y, z) — (72, 0, 0) with 7 the constant shear 
rate. 

The initial state was taken to be a sub-monolayer of 
particles with number density per unit surface area. 
Particles were added at random uniformly over the sur- 
face, and attempted additions that would create particles 
with overlapping cell radii were rejected. Each seed parti- 
cle was anchored to a point directly beneath it (see below 
for the definition of anchors). 

Although the three model components share the same 
spatial domain, they relax on separated time scales, al- 
lowing them to be solved sequentially. The iteration 
cycle proceeds in the order scalar — hiomass — fluid 
— >■ scalar with data extraction just before the 

hiomass growth iteration. Each stage in this cycle is now 
explained in detail. 




(a) (b) 

FIG. 3: (a) Single particle i with cell diameter d1 and EPS 
diameter dl. (b) Schematic of redistribution of cellular and 
EPS masses after particle i (dashed lines) divides into i\ and 
12. Each mass component is conserved. 



B. Scalar iteration 

The nutrient concentration c(x) obeys the steady-state 
reaction-diffusion-advection equation 

N 

= 5tC = -v-Vc + L>V2c-t-^r,(5(x-x,) (2) 

i=l 

obeying the mixed boundary conditions 

c{z = L^) ^ Co, (3) 
a.cL^o - 0. (4) 

The reaction rates are given by Q, and note that dif- 
fusion is assumed to be constant. This is solved numeri- 
cally using a finite difference scheme solved on a regular 
rectangular mesh using geometric multigrid |29j . To de- 
termine the reaction terms in ([2]), the value of c at the 
particle centre x^ is found by trilinear interpolation from 



Label 


Meaning 


Value 


Co 


Bulk nutrient concentration 


- 


7 


Fluid shear rate 




Lx 


Box length in direction of flow 


- 


Ly 


Box width in vorticity direction 






Height from solid surface to bulk 


400 Atm 




Cell density (excluding water) 
EPS density (excluding water) 


0.2 pg/^m^ 
4x10"^ Pg/MHi^ 


Ki/2 


Monod half-concentration 


10"® Pg/^trn^ 
10^ ^mVs 


D 


Nutrient diffusion coefficient 


^max 


Base reaction rate 


0.5/h 


yc 


Yield factor for cell mass 


0.2 


rcl 


Relative yield factor for EPS 


0.4 


^max 


Division diameter 


5^m 


div 

(J 


Width of relative mass division 


n 1 
U.l 


V 


Fluid viscosity 


10"^ Pa s 


^anc 


Anchor spring stiffness 


50 pN/^m 




EPS spring stiffness per mass 


5 pN [ivcT^ Pg~^ 




Initial surface number density 


10"^ /im"^ 



TABLE I: Variables and parameters. Those above the line are 
treated as variables here, while those below were kept fixed 
with the values quoted, which were chosen to be representa- 
tive of oral bacteria taking simple sugars as a nutrient |26H28) . 



the adjacent mesh nodes, inserting into ([T]), and then dis- 
tributing the resulting onto the same lattice nodes in 
a way that conserves the total reaction rate. Here we 
weight the contribution to each node by the inverse of its 
distance from x^. 



C. Biomass iteration 

Once the steady-state reaction rates have been de- 
termined for each particle i, the increase in both cellu- 
lar mass ml and the EPS mass m\ are found by mul- 
tiplying the mass growth rates by the biomass time in- 
terval At'"". This time-step is adaptive, so that higher 
relative growth rates correspond to smaller time steps 
and vice versa. A linear variation was employed here, 

At^'° = Cmaxi=i,„N (^^r^tmi), with C = 0.01 to fix 

the maximum particle growth at around 1% per time 
step. C was varied to ensure no discernible variation 
of measured quantities. Note that this biomass growth 
time step is distinct from, and many orders of magnitude 
larger than, the time step used during fluid stabilisation 
described below. 

After the cellular and EPS masses of each particle, 
and thus their corresponding diameters, have been up- 
dated, the system is checked for division events. Any 
particles whose new diameter exceeds the division thresh- 
old, i.e. dl > d'^^^, divides into two daughter cells ii 
and i2. Mass is conserved during division, but is dis- 
tributed asymmetrically between the two daughters ac- 
cording to TO^^ = m1 — m1^ = XiTTil, where is a random 
variable chosen for each division event from a Normal dis- 
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tribution with mean | and width cr'^'^. The EPS mass is 
divided similarly, with the same A^. The daughter cells 
are placed at opposite poles of a sphere, centred on the 
parent cell, with a diameter D = ^ (d^^ ) + ^ (dj^ ) 
so that their EPS shells overlap; see Fig.jSTb). The axis of 
the sphere on which the daughters are added is chosen at 
random to ensure division cannot introduce anisotropy. 

The links between the particles can now be determined. 
In essence, this amounts to identifying pairs of particles 
i and j whose EPS shells overlap, |xi — Xj| < ^(df + d'j), 
and adding a spring between their centres. In practice 
this leads to the rare instances where both daughter cells 
become disconnected from the film shortly after a di- 
vision event. To robustly maintain film integrity, after 
each round of division events, all particles are sorted into 
clusters, where two particles belong to the same clus- 
ter if their EPS shells overlap. Any isolated clusters are 
translated into contact with either the film or the base at 
z = 0, whichever requires the shortest motion. Further- 
more, to avoid 'knots' of springs, no particle is allowed 
to have more than 13 links. Any particle with more than 
this number of links has its longest links removed until 
this maximum number is reached. 

Links between particle pairs are deleted before each 
growth and division cycle, and recreated afterwards. 
They are therefore transient links that reflect the cur- 
rent configuration of the film. Links to the base at z = 
are different in that they cannot move once formed, else 
the film could drift in an uncontrolled manner in the pres- 
ence of flow. Instead, these anchor links are permanent 
and do not move once formed. They are created when a 
particle that does not already have an anchor link comes 
into contact with the base, i.e. has a centre at a height 
Zi < ^d°. A spring is then created between the particle 
and an anchor point that is directly below the particle 
position at this time, i.e. at (xi, 2/^,0). An anchor is 
not created if the particle already has 3 transient links to 
anchored particles. These rules maintain a stable popula- 
tion of anchor links that does not drift during the biofilm 
evolution. 



D. Fluid iteration 

In a full model with a spatio-temporally varying flow 
field, v(x) would need to be simultaneously solved 
with the stabilisation of the biomass in a momentum- 
conserving manner. Since v(x) is fixed here, we need only 
stabilise the film in the presence of flow. This amounts to 
demanding that the net force on each particle i simul- 
taneously vanishes. Two forces contribute to fi, a drag 
force deriving from the flow, and a matrix force due to 
the links between particles, or anchor links to the base. 
The drag force is based on Stokes flow past a sphere, 

ff^s = 3W?v(x,) (5) 

where the fluid viscosity v is chosen to be that of wa- 
ter. The matrix force f™*** derives from the links de- 



termined in the previous step that are now identifled as 
Hookean springs. For anchor links, the spring force is 
^ancj-j, — id^) where r is the distance of the cell cen- 
tre from the anchor point on the surface, and fc^""^ is 
a uniform spring constant. This scalar force is pro- 
jected along the line connecting the particle to the anchor 
point to give the required vector force. For the tran- 
sient EPS-mediated links between particle pairs i and j, 
the force is K°TO°j(|xi — Xj \ — £q) with a natural length 
io = \{di + d'j + d^ + d'j) corresponding to the midpoint 
of the EPS shells. Here, k" is the stiffness per unit mass, 
and m^, is the mass of the EPS that is attributed to this 
link. This is determined by equally distributing each par- 
ticle's EPS mass to each of its (non-anchor) links. This 
scalar force is projected to the line of centres between x^ 
and Xj in an equal-and-opposite manner. 

The goal is to determine the whole film configuration 
{xj for which each f, = ff'^'s + fmat ^ q_ rp^^ 
ods were used here which gave equivalent results. Since 
they are standard they will only be described briefly here. 
The non-linear conjugate gradient method [SIT, which was 
found to be most efficient for small systems, requires 
repeated construction and inversion of a large, sparse 
stiffness matrix giving the changes in each component of 
each force for small changes in particle positions. Block- 
diagonal preconditioning was also used. The second 
method, which proved to be more efficient for large sys- 
tems and those with fiow, was to use overdamped molec- 
ular dynamics [BTj in which particles were moved in the 
direction of their unbalanced force: Ax^ = Atfi/Siri/d'i, 
where an adaptive time step At was used that increases 
as the largest velocity decreases. For both methods, con- 
vergence tolerances were systematically varied until there 
was no discernible variation in measured quantities. 



III. RESULTS 

The control variables are here chosen to be those that 
are also amenable to experimental control, namely the 
bulk nutrient diffusion Cq and the shear rate 7. The hor- 
izontal surface dimensions = Ly are systematically 
varied to determine finite size effects. All other param- 
eters are kept fixed with the values quoted in Table [Tj 
which were taken to be representative of oral bacteria 
growing in the presence of sugars p6H28] . The theoretical 
predictions for fiat films referred to below are derived in 
Appendix]^ Unless otherwise stated, all results are pre- 
sented in terms of dimensionless quantities constructed 
by scaling by combinations of the length d"^^^, inverse 
time /Cmax and mass concentration ^^1/2- 



A. Surface roughening without flow 

We start with the no-fiow case 7 = 0.0. The mean 
surface height h{t) was extracted from the simulations 
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using the procedure explained in Appendix [B] For all pa- 
rameters studied, the variation of h{t) with time showed 
no significant variation with the horizontal system size 
Lx = Ly. An example is given in Fig.|4j where the analyt 



ical solution for a flat film ( Af 1 is also plotted. The bulk 



cell mass density nm in this solution curve was measured 
independently, so there are no fitting parameters. Actual 
growth curves consistently exceed this theoretical predic- 
tion at late times. Since the degree of overshoot corre- 
lates with increasing surface roughness (as defined below; 
data not given), we infer this results from the omission of 
surface undulations in the calculations. Note that direct 
observation of the data confirms that c{z < h) <^ K1/2 
in all cases, as per the calculations. 





FIG. 5: Surface roughness versus height for co = IQK1/2, 
7 = and the horizontal system sizes given in the legend. 
For comparison, the solid black line segment has a slope of 
~ 0.27. Bars show standard error over independent runs with 
differently randomised initial conditions and mass redistribu- 



tion after division (n = 
30d'^=''=, AQcT'^'' resp.). 



-18, 10, 6 runs for Lx 



20d" 



The influence of finite system size is expected to be 
due to a horizontal correlation length that grows with 
time, causing w to saturate when approaches ~ Ly. 

can be extracted from the height-height correlation 
function Chhir), 



FIG. 4: Mean surface height h versus time t for cq = 5/4'i/2 
and no flow, 7 = 0. The solid black line shows the flat film 



prediction from (AlOl, which requires no fitting parameters. 



The height has been scaled to the threshold diameter for par- 
ticle division d™^", and t to the base reaction rate fcmax. The 
horizontal system sizes are given in the legend. The series for 
the largest system is shorter due to computational limitations. 

As the mean height h{t) grows non-linearly with time, 
unlike the canonical surface growth models where it 
grows at a constant rate [TT], we hereafter take ft, as a 
surrogate time variable to permit direct comparison with 
other models. We first consider the surface roughness or 
width w defined by 



w 



LxLy 



dx / dy [h{x,y) - h]' 



(6) 



where h{x, y) is the height of the surface vertically above 
basal coordinates {x, y). The typical growth of w with h 
is shown in Fig. [Sj As with fractal growth models, this 
increases until saturating at a maximum value that in- 
creases with system size. Unlike canonical models, where 
growth is sub- linear |TT], here the growth rate is consis- 
tent with Zmear scaling w(t) cx h{t) as shown in the figure. 
Unfortunately the poor statistics rules out a more precise 
evaluation of the growth exponent. 



Chkir) = dx dy [h{0, 0)h{x, y) ~ 



(7) 



where r'^ — x'^ + y^ and translational symmetry in the x- 
y plane has been assumed. For all plots, Chhir) crossed 
from positive at small r to negative at large r (data not 
shown); the single crossing point is identified with . An 
exaniple of the variation of ^" and system size is given in 
Fig. |6j and confirms the expected picture that ^" grows 
with time until reaching a maximum value that increases 
with system size. The variation of with h before satu- 
ration is again consistent with linear growth as shown in 
the figure, and the data for other parameters, although 
noisier, are consistent with this growth law. 

In [23] it was observed that rougher films correlated 
with a thinner layer of actively growing particles near the 
surface, as compared to a thick active layer that gener- 
ated flatter films. In [23 the thickness of the active layer 
was quantified by an a priori function of input parame- 
ters. Here we instead directly measure the layer thick- 
ness, and compare to the measured roughness w. The 
active layer is defined in terms of the relative growth rate 
of particles {m'i)~^dtm1 measured as a function of verti- 
cal distance Az from the surface. Details of how this was 
extracted from the simulations is given in Appendix [B] 
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FIG. 6: Growth of horizontal correlation length for the system 
sizes given in the legend and the same parameters as Fig. [S] 
The straight black line segment has a slope of 0.24. 



FIG. 7: Depth of the active layer defined by Q versus height 
for the same parameters as Fig. [S] 



The penetration depth is then 



I Az=0 



(a.(K)-i5,m?))|, 



(8) 



where the averaging (• ■ ■ ) is over all particles at the 
same depth Az > below the surface. The variation 
of is plotted in ([T]) and demonstrates weak growth. 
By contrast, the flat-film prediction (A6) takes a value 
« 1.85 c?™'''^, 20-30% smaller than measured, and does 
not increase with time. Again, the likely culprit for the 
excess measured thickness is the inapplicability of the 
flat-film assumption. Note that although the theoretical 
prediction employs the variation of c(z) rather than the 
growth rate, c is roughly proportional to growth for the 
considered parameters, so these two definitions of 1^ are 
equivalent. 

To correlate with roughness, it is convenient to re- 
duce both time-varying quantities to single scalars. For 
the roughness, we focus on the linear growth regime 
w = ah + b and extract the slope a as a measure of 
roughness. For the active layer, we take the average of 
£p over the region of slow growth shown in Fig. [Tj in 
the understanding this is just a working definition that 
will weakly depend on the achievable simulation times. 
Plotting these two as in Fig. [8] shows an inverse relation- 
ship between roughness and the depth of the active layer, 
confirming the finding of [23] . 



B. Effect of afflne flow 

We now turn to consider the effects of flow, 7 > 0, 
keeping the bulk concentration fixed at cg = 10^^1/2- Al- 
though the mean surface height grows at a slightly lower 
rate in the presence of flow, much more striking is the 




FIG. 8: Surface roughness versus depth of the active layer. 
Closed symbols correspond to 7 = 0, with increasing Co/ K\/2 
— 1, 5, 10 and 20 as indicated by the upper arrow. Open 
symbols correspond to cq = 107^1/2 and increasing 7 = 0, 
0.072 fcniax and 0.72 fcmax as indicated by the lower arrow (the 
7 = point belongs to the connecting point in the first data 
set). 



significant decrease in surface roughness demonstrated in 
Fig. [9] Although the growth law remains approximately 
linear, the slope is noticeably reduced compared to the 
no-flow case. It might be postulated that the reduction 
in roughness is due to some change in the depth of the 
active layer. However, as shown in Fig.|8] flow affects the 
roughness but not the depth of the surface layer. 

Systematically varying the system size reveals a mixed 
picture. For the highest flow rate 7 « 0.72fcniax, there is 



no significant variation with = Ly as shown in Fig. 
For the lower flow rate considered, 7 « 0.072 fcmax, 
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roughness for = Ly = 20 d™'^'^ significantly exceeded 
that for Lx = Ly = 30 d™^^, taking values close to the 
7 = case, although the statistics are poor and this ob- 
servation is not definitive. This uncertainty is reflected in 
the large vertical error bar for this point in Fig. [8] While 
the reduction of roughness due to shear is clear, improved 
statistics and a larger range of system sizes, both requir- 
ing the development of more advanced algorithms, will 
be required to fully clarify the picture. 

A final observation relates to the mean cellular mass 
density, denoted nm in connection with the theory of 
Appendix [Xj This was measured for all parameters and 
system sizes far from the surface, and exhibited no sig- 
nificant variation with system size or cq. It did however 
admit a slight but definite decrease for high flow rates, 
dropping roughly 5% for the highest flow rate considered, 
7 w 0.72fcinax, compared to 7 = 0. This is most prob- 
ably an expression of Reynolds' dilation, a phenomenon 
common to particulate media where shear stresses gener- 
ate system-spanning force chains that react against the 
solid surface, raising the system and lowering the mean 
density [32j . 













^fn 7s:0.072 fc„,„ 

7-0.72 fc„,„ 



h/d" 



FIG. 9: Growth of roughness for the shear rates given in 
the legend, for system size Lx = Ly — 20 d"^^^. Note that 
the data sets are shorter for the fastest flow rate considered 
here as the mechanical stabihsation algorithm staUed for thick 
films, necessitating premature termination of the simulation. 



IV. DISCUSSION 

Many features common to fractal growth models [TT] 
have been observed in this investigation, including an 
algebraic increase in surface roughness in both the hor- 
izontal and vertical directions, that saturates when the 
horizontal correlations ^" become comparable to the sys- 
tem size. Somewhat anomalous are the growth expo- 
nents themselves, which are both consistent with linear 
growth, significantly faster than the sub-linear laws typ- 
ically measured. An explanation based on 'freezing' of 
surface regions might provide a simple explanation for 
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FIG. 10: Variation of the growth in surface roughness with 
system size for 7 ~ 0.072 fcmax. Compare to the no-flow case 
in Fig. [5] 



the linear growth in w, but not in ^H, and in any case is 
not consistent with direct observation of the full surface 
profiles which suggest no such freezing. Although the 
linear growth of w and suggests a dynamic exponent 
also equal to I [11] , our statistics are too poor to permit a 
meaningful check of this additional exponent to confirm 
this relation. 

We hypothesise that a key mechanism underlying the 
rapid increase in roughness is the non-locality deriving 
from the nutrient concentration field c(x). In this con- 
text, it is instructive to note that the stationary Green's 
function (i.e. the steady solution for a point source) 
for ll), in an infinite system without flow, decays with 
distance |x| from the source as c(x) c>c |x|~^ [33], a scale- 
invariant, long-range decay. In 2D the same solution does 
not decay at all but increases logarithmically. Non-local 
effects should therefore be expected. Nonetheless we have 
been unable to derive a simple explanation for the w oc h 
growth law, and suggest that the construction of simpli- 
fied models will allow larger systems to be reached and 
generate insight into this phenomenon. Furthermore, we 
cannot rule out a crossover to different scaling at late 
times exceeding our simulation capabilities, as in some 
other models with non-local surface interactions [TTJ [T^] . 
We note however that the biofilm thickness reached in 
our simulations, roughly 150 — 200/im, are comparable to 
real biofilms, therefore our findings should be regarded as 
biologically relevant. 

Shear flow is well known to induce waves at liquid sur- 
faces, but can also smooth surfaces by suppressing ther- 
mal capillary waves as observed in colloidal gas-liquid 
interfaces [34|. This insight cannot be transferred to our 
athermal system, however, thus the mechanism under- 
lying the measured reduction in roughness is not clear. 
Since the elastic strains in the biofilm were visibly very 
small, the observed smoothing is most likely due to alter- 
ations to the transport of the nutrient. It is not simply 
due to changes to the mean nutrient transported to the 
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surface, however, as this would affect the depth of the 
active layer, which was not observed. Instead, we argue 
that the effect of flow on roughness can be intuitively 
understood as being due to the competition between dif- 
fusion, controlled by the parameter D, and advection due 
to the velocity field v(x) = (72,0,0). The effect of dif- 
fusion over advection can be quantified by the dimen- 
sionless Prandtl number P ~ D/{jh'^) with h a char- 
acteristic height of the film. Taking h w lOOfim as a 
typical biofilm thickness, the two values of 7 employed 
here, 7 « O.G72fcinax and 0.72fcinax, correspond to P « 10 
and P « 1 respectively. This confirms the relevant role of 
advection to our results. It does not, however, highlight 
the microscopic mechanism underlying the smoothing, 
and here again further investigation of simplified models 
is desirable. 

This first application of the mechanical-IbM model re- 
mains deficient in two key respects. Firstly, the coupling 
between the fluid and the rest of the system is strictly 
one-way, i.e. the fluid affects the biofilm and the nutrient 
field, but the biofilm docs not affect the flow. Intuitively 
we would expect the fluid to flow around surface features, 
potentially unlocking a rich interplay between nutrient 
transport and surface morphology. Secondly, biomass 
detachment due to shear stresses has not been incorpo- 
rated, although this is known to partly control biofilm 
thickness [351 136j . There is no reason why this cannot 
be introduced for a future work, possibly following the 
particle-removal criterion of Alpkvist and Klapper [37]. 
We note that the m-IbM model maintains the primary 
advantages of the IbM approach, including the relative 
ease of modelling multi-species films, and speculate it 
will become an important tool in the quantification of 
biofilm-flow coupling in the future. 
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Appendix A: Flat-film theory 

For flat films with uniform thickness h{t), it is pos- 
sible to write down analytically tractable equations by 
approximating the biofilm as a continuous body. The 
discrete reactions in ([l]) are replaced by the continuous 
field r{z), which is proportional to the number density of 
cells per unit volume, n, and the mass per cell, m, both 
of which are taken as uniform and constant. For clar- 
ity of the resulting expressions, we define a = nmk^^^. 
Then c(z) obeys the following one-dimensional reaction- 



diffusion equation, 

{)^dtc{z) = Ddlc{z) + r{z), 


c{z) 



r{z) 



-a 



C{z) + 



z > h{t), 
z < h{t). 



(Al) 
(A2) 



Even with these simplifications, (A2| is non-linear and 
no general analytical solution is apparent. Instead we 
consider limits of high and low c throughout the biofilm, 
i.e. c{z < h) ^ Ki/2, and conversely c{z < h) < C K1/2, 
for which the non-linearity is removed and (Al) can be 
solved. The solution for c{z < /i) ^ K1/2 is 



co-c(z) = 



D 

ah 



P. - 



2D 



> h. 



< h. 



(A3) 



For consistency we must also have c(0) ^ K1/2. In the 
opposite limit c(z < /i) <C K1/2, and for clarity defining 

13^ = a{Ky2Dr\ 



c{zl 
Co 



1 - 



/3(P^ - z) smh{ph) 



cosh(/3/i) + P{L, - h) sinh(/3/i) 

cosh(/3z) 
^ cosh(^/i) -I- /3(P^ - h) sm\i{j3h) 



:z > h, 



z < h. 



(A4) 

Here we additionally require c{h) <^ K1/2. For both 
limits, continuity of c(z) and dzc{z) at z — h{t) has been 
imposed. 

For c{z < /i) 3> Ki/2 > the concentration remains 
significantly non-zero to the base of the film. By con- 
trast, for c{z < h) <^ Ki/2 the concentration can become 
vanishingly small while still within the film. In this case 
we define the penetration depth by 



a.c(z)L 



(A5) 



where continuity of the first derivative means that either 
of the z < h or z > h expressions in (A4) can be used, 
giving 



coth(/3/i) 



(3-^ for > 1. 



(A6) 



For thick films /3h ^ 1 this increases with D and de- 
creases for higher reaction rates. It is thus a length scale 
that determines the balance of diffusion to reaction, and 
plays a comparable role to the (dimensionless) Thiele 
modulus [38J. Conversely for thin films /3ft- <C 1, di- 
verges, suggesting it cannot be identified with a physical 
length scale in the original discrete system. 

The time evolution of the film thickness h{t) can be 
determined by considering the rate of change of the to- 
tal mass in the film, and maintaining the assumption of 
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constant nm. It is then straightforward to derive the fol 
lowing integro-differential equation from (All and (A2) 



dh{t) 
At 



h{t) 



dz 



c[z) 



C{z) + 



(A7) 



It is again simpler to remove the non-linearity in the inte- 
grand by considering limits of c(z). For c(z < h) ^ Kij2, 
when the nutrient penetrates throughout the entire film, 
(A7) is readily solved to give exponential growth, 



hit) = h{Q)t 



(A8) 



Finally, note that the crossover between shallow and deep 
penetration can be expressed in terms of the dimension- 
less ratio ^p/ (L^ — h), which (for ip « /S^^) gives a similar 
quantity to the S employed in |23) . 



Appendix B: Data analysis of the surface 

To determine the moments of the height distribution it 
is first necessary to identify the surface. To do this, the 
system box was partitioned into a cubic lattice in which 



For c(z < h)^ K,/2 (|A4j) can be used to evaluate the ^^^^ ^^^^^ dimensions d 



integral in ( A7 1 , producing the differential equation 



Every 



dh{t) _ coDY" 
dt ^ 



1 



nm ip + (L^ - h{t)) 



Note that there is implicit /i-dependence in the nutri- 
ent penetration depth ip as per (A6). Because of this. 



it is simplest to solve (A9) in the limits of shallow and 



deep nutrient penetration layers relative to the boundary 



layer. 



<^ Lz — h and £p ^ L^ — h respectively. For 



the former case, the solution is 



h{t) = L, 



2coDY^ 



t 



This expression predicts the film reaches Lz at a fi- 
nite time [Lz — h{0)]'^nm/ {2cqDY'^) , but the assumption 
Ip ^ Lz — h{t) will break-down before this happens. The 
corresponding solution for the deep penetration depth 
h is 



limit £p > L 



h{t) = /3"^arsinh \ smh[l3h{0)] exp 



V ^1/2 



lattice block with one or more particle centres con- 
tained within it was marked as occupied; all others are 
(A9) marked vacant. Lattice blocks on the base, i.e. in the 
plane z = 0, are labelled k and I in the x and y-directions 
respectively. The height hki of the film above each base 
block is defined as the midpoint of the highest occupied 
block vertically above it. Note that this definition ig- 
nores overhangs, but as these were rarely observed they 
should only represent a small correction to our basic find- 
ings. Moments of hjk were calculated as per any discrete 
distribution. For the spatial correlations in height, the 
horizontal distance r between midpoints of base lattice 
(AlO) blocks were used, incorporating the periodic boundary 
conditions in the horizontal directions. 

Metabolic activity as a function of distance from the 
surface was measured using the same lattice. In this case, 
the mean relative growth rate •m~^dtm of all particles in 
each lattice block were calculated and assigned to that 
block. This was output as a function of distance from 
the highest occupied site in the same column (fc,?), so 
a depth of corresponds to the growth rate of particles 
in the highest occupied block, rf™^'^ to the block directly 
(All) beneath it, etc. 
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